##rm(list = ls()) 


install.packages("foreign") 
install.packages("data.table")
install.packages("car")
install.packages("readstata13")
install.packages("plotrix")
install.packages("descr")
install.packages("weights")
install.packages("miceadds") 
install.packages("mice")
install.packages("readxl")
install.packages("scales")
install.packages("ggplot2")
install.packages("gridExtra")
install.packages("grid")
install.packages("MASS")
install.packages("Hmisc")
install.packages("qte")
install.packages("cem")
install.packages("MatchIt")
install.packages("Zelig")
install.packages("RItools")
install.packages("dplyr")
install.packages("poliscidata")
install.packages("tidyr")
install.packages("sandwich")
install.packages("lmtest")

library(foreign) 
library(data.table)
library(car)
library(readstata13)
library(plotrix)
library(descr)
library(weights)
library(miceadds) 
library(mice)
library(readxl)
library(scales)
library(ggplot2)
library(gridExtra)
library(grid)
library(MASS)
library(Hmisc)
library(qte)
library(cem)
library(MatchIt)
library(Zelig)
library(RItools)
library(dplyr)
library(poliscidata)
library(lmtest)
library(tidyr)
library(sandwich)
library(lmtest)

# # Set the working directory to the location of the current script
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))  # Set working directory to the current script's directory
getwd()  # Check the current working directory

# load in file when it isnt in the environment already
load("HurricaneDataFinal.RData")



#================================#
# TREATMENT COUNTY RELIEF POLICY #
#================================#

# Voted2016 vs. Voted 2018 (any method, early or election day)

# Reshape dataset from wide to long
dta_ReliefSurrounding_long <- reshape(dta_ReliefSurrounding,
                                      varying = c("Voted2016", "Voted2018"),
                                      v.names = "voted",
                                      timevar = "year",
                                      times = c("Voted2016", "Voted2018"),
                                      new.row.names = 1:946286,    #multiply the number of obs by 2
                                      direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurrounding_long$year2018 <- ifelse(dta_ReliefSurrounding_long$year=='Voted2018', 2018, 2016)
dta_ReliefSurrounding_long$d18 <- ifelse(dta_ReliefSurrounding_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurrounding_long <- data.table(dta_ReliefSurrounding_long)


# DiD variable: treatment*year
dta_ReliefSurrounding_long$did <- dta_ReliefSurrounding_long$Relief * dta_ReliefSurrounding_long$d18

#table(dta_ReliefSurrounding_long$did)



# Difference-in-Difference County Treatment Model Estimation 
# m_ReliefSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurrounding_long)
# 
# summary(m_ReliefSurrounding)

m_ReliefSurrounding <- coeftest(m_ReliefSurrounding, vcov=vcovHC(m_ReliefSurrounding, type="HC0", cluster="CountyCode.x")) 

m_ReliefSurrounding

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurrounding_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.9674 -0.4628  0.1628  0.3377  0.8609 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6491697  0.0015084  430.364  < 2e-16 ***
#   d18         -0.1122619  0.0012302  -91.257  < 2e-16 ***
#   Relief      -0.0567112  0.0012988  -43.664  < 2e-16 ***
#   did         -0.0722536  0.0018273  -39.541  < 2e-16 ***
#   AgeCat20182  0.0933976  0.0013996   66.731  < 2e-16 ***
#   AgeCat20183  0.2560783  0.0013259  193.129  < 2e-16 ***
#   AgeCat20184  0.3098964  0.0014351  215.939  < 2e-16 ***
#   gender2     -0.0235947  0.0009196  -25.656  < 2e-16 ***
#   Dem         -0.0473426  0.0011234  -42.144  < 2e-16 ***
#   NPA         -0.1947494  0.0013538 -143.850  < 2e-16 ***
#   Third       -0.0744080  0.0064010  -11.624  < 2e-16 ***
#   Black        0.0083583  0.0013076    6.392 1.64e-10 ***
#   Hispanic    -0.0207154  0.0028978   -7.149 8.77e-13 ***
#   Other       -0.0505346  0.0021415  -23.598  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.441 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.127,	Adjusted R-squared:  0.127 
# F-statistic: 1.052e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# i.e. those living in counties hit by michael who voted in 2016 (by any method) were less likely to vote (by any method) in 2018 than were those living in counties not hit by michael who voted in 2016 (by any method)


### VOTING METHODS ###


# Voted 2016 vs. Any Early Method 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted16VsAnyEarly_long <- reshape(dta_ReliefSurrounding,
                                                       varying = c("Voted2016", "AnyEarlyMethod2018"),
                                                       v.names = "voted",
                                                       timevar = "year",
                                                       times = c("Voted2016", "AnyEarlyMethod2018"),
                                                       new.row.names = 1:946286,    #multiply the number of obs by 2
                                                       direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVoted16VsAnyEarly_long$year2018 <- ifelse(dta_ReliefSurroundingVoted16VsAnyEarly_long$year=='AnyEarlyMethod2018', 2018, 2016)
dta_ReliefSurroundingVoted16VsAnyEarly_long$d18 <- ifelse(dta_ReliefSurroundingVoted16VsAnyEarly_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted16VsAnyEarly_long <- data.table(dta_ReliefSurroundingVoted16VsAnyEarly_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted16VsAnyEarly_long$did <- dta_ReliefSurroundingVoted16VsAnyEarly_long$Relief * dta_ReliefSurroundingVoted16VsAnyEarly_long$d18

table(dta_ReliefSurroundingVoted16VsAnyEarly_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted16VsAnyEarlySurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted16VsAnyEarly_long)
# 
# summary(m_ReliefVoted16VsAnyEarlySurrounding)

m_ReliefVoted16VsAnyEarlySurrounding <- coeftest(m_ReliefVoted16VsAnyEarlySurrounding, vcov=vcovHC(m_ReliefVoted16VsAnyEarlySurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted16VsAnyEarly_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.9734 -0.4198  0.1096  0.3854  0.9476 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6593099  0.0015131  435.727  < 2e-16 ***
#   d18         -0.3730352  0.0012340 -302.294  < 2e-16 ***
#   Relief      -0.0543818  0.0013029  -41.740  < 2e-16 ***
#   did          0.0429501  0.0018330   23.431  < 2e-16 ***
#   AgeCat20182  0.0666737  0.0014040   47.489  < 2e-16 ***
#   AgeCat20183  0.2202640  0.0013301  165.602  < 2e-16 ***
#   AgeCat20184  0.3087043  0.0014396  214.439  < 2e-16 ***
#   gender2     -0.0232394  0.0009225  -25.192  < 2e-16 ***
#   Dem         -0.0446825  0.0011269  -39.652  < 2e-16 ***
#   NPA         -0.1637920  0.0013581 -120.607  < 2e-16 ***
#   Third       -0.0649231  0.0064210  -10.111  < 2e-16 ***
#   Black        0.0053499  0.0013117    4.079 4.53e-05 ***
#   Hispanic    -0.0081587  0.0029068   -2.807    0.005 ** 
#   Other       -0.0354116  0.0021482  -16.485  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4424 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.2055,	Adjusted R-squared:  0.2055 
# F-statistic: 1.871e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# i.e. those living in counties hit by michael who voted (by any method) in 2016 were MORE likely to vote early (by any method) in 2018 than were those living in counties not hit by michael who voted (by any method) in 2016 


# Voted 2016 vs. Election Day 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted16VsED_long <- reshape(dta_ReliefSurrounding,
                                                 varying = c("Voted2016", "EDOnly2018"),
                                                 v.names = "voted",
                                                 timevar = "year",
                                                 times = c("Voted2016", "EDOnly2018"),
                                                 new.row.names = 1:946286,    #multiply the number of obs by 2
                                                 direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVoted16VsED_long$year2018 <- ifelse(dta_ReliefSurroundingVoted16VsED_long$year=='EDOnly2018', 2018, 2016)
dta_ReliefSurroundingVoted16VsED_long$d18 <- ifelse(dta_ReliefSurroundingVoted16VsED_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted16VsED_long <- data.table(dta_ReliefSurroundingVoted16VsED_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted16VsED_long$did <- dta_ReliefSurroundingVoted16VsED_long$Relief * dta_ReliefSurroundingVoted16VsED_long$d18

table(dta_ReliefSurroundingVoted16VsED_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted16VsEDSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted16VsED_long)

# summary(m_ReliefVoted16VsEDSurrounding)

m_ReliefVoted16VsEDSurrounding <- coeftest(m_ReliefVoted16VsEDSurrounding, vcov=vcovHC(m_ReliefVoted16VsEDSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted16VsED_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.85065 -0.24089 -0.05102  0.24841  1.09366 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.7073516  0.0014168  499.262  < 2e-16 ***
#   d18         -0.4871646  0.0011555 -421.623  < 2e-16 ***
#   Relief      -0.0405730  0.0012199  -33.259  < 2e-16 ***
#   did         -0.0922082  0.0017163  -53.724  < 2e-16 ***
#   AgeCat20182  0.0640476  0.0013146   48.720  < 2e-16 ***
#   AgeCat20183  0.1433016  0.0012454  115.064  < 2e-16 ***
#   AgeCat20184  0.1313528  0.0013479   97.447  < 2e-16 ***
#   gender2     -0.0174272  0.0008638  -20.176  < 2e-16 ***
#   Dem         -0.0269462  0.0010551  -25.539  < 2e-16 ***
#   NPA         -0.1225683  0.0012716  -96.389  < 2e-16 ***
#   Third       -0.0470251  0.0060122   -7.822 5.22e-15 ***
#   Black       -0.0094365  0.0012282   -7.683 1.55e-14 ***
#   Hispanic    -0.0163964  0.0027218   -6.024 1.70e-09 ***
#   Other       -0.0410678  0.0020114  -20.418  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4142 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.3117,	Adjusted R-squared:  0.3117 
# F-statistic: 3.275e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted (by any method) in 2016 were less likely to vote on Election Day in 2018 than were those living in counties not hit by michael who voted (by any method) in 2016

# Voted 2016 vs. VBM Only 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted16VsVBM_long <- reshape(dta_ReliefSurrounding,
                                                  varying = c("Voted2016", "VBMOnly2018"),
                                                  v.names = "voted",
                                                  timevar = "year",
                                                  times = c("Voted2016", "VBMOnly2018"),
                                                  new.row.names = 1:946286,    #multiply the number of obs by 2
                                                  direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVoted16VsVBM_long$year2018 <- ifelse(dta_ReliefSurroundingVoted16VsVBM_long$year=='VBMOnly2018', 2018, 2016)
dta_ReliefSurroundingVoted16VsVBM_long$d18 <- ifelse(dta_ReliefSurroundingVoted16VsVBM_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted16VsVBM_long <- data.table(dta_ReliefSurroundingVoted16VsVBM_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted16VsVBM_long$did <- dta_ReliefSurroundingVoted16VsVBM_long$Relief * dta_ReliefSurroundingVoted16VsVBM_long$d18

table(dta_ReliefSurroundingVoted16VsVBM_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted16VsVBMSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted16VsVBM_long)

# summary(m_ReliefVoted16VsVBMSurrounding)

m_ReliefVoted16VsVBMSurrounding <- coeftest(m_ReliefVoted16VsVBMSurrounding, vcov=vcovHC(m_ReliefVoted16VsVBMSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted16VsVBM_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.91535 -0.19264 -0.01512  0.23714  1.12843 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6931649  0.0012839  539.884  < 2e-16 ***
#   d18         -0.6142461  0.0010471 -586.626  < 2e-16 ***
#   Relief      -0.0446918  0.0011055  -40.426  < 2e-16 ***
#   did         -0.0014815  0.0015554   -0.953   0.3408    
# AgeCat20182  0.0452023  0.0011913   37.944  < 2e-16 ***
#   AgeCat20183  0.1383371  0.0011286  122.574  < 2e-16 ***
#   AgeCat20184  0.2221850  0.0012215  181.893  < 2e-16 ***
#   gender2     -0.0246202  0.0007828  -31.453  < 2e-16 ***
#   Dem         -0.0239508  0.0009562  -25.049  < 2e-16 ***
#   NPA         -0.1094763  0.0011523  -95.003  < 2e-16 ***
#   Third       -0.0387170  0.0054483   -7.106 1.19e-12 ***
#   Black       -0.0142527  0.0011130  -12.805  < 2e-16 ***
#   Hispanic    -0.0058276  0.0024665   -2.363   0.0181 *  
#   Other       -0.0270825  0.0018227  -14.858  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3754 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.4252,	Adjusted R-squared:  0.4252 
# F-statistic: 5.349e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted (by any method) in 2016 were MORE likely to vote by mail in 2018 than were those living in counties not hit by michael who voted (by any method) in 2016

# Voted 2016 vs. EIP Only 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVoted2016VsEIP_long <- reshape(dta_ReliefSurrounding,
                                                    varying = c("Voted2016", "EIPOnly2018"),
                                                    v.names = "voted",
                                                    timevar = "year",
                                                    times = c("Voted2016", "EIPOnly2018"),
                                                    new.row.names = 1:946286,    #multiply the number of obs by 2
                                                    direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVoted2016VsEIP_long$year2018 <- ifelse(dta_ReliefSurroundingVoted2016VsEIP_long$year=='EIPOnly2018', 2018, 2016)
dta_ReliefSurroundingVoted2016VsEIP_long$d18 <- ifelse(dta_ReliefSurroundingVoted2016VsEIP_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVoted2016VsEIP_long <- data.table(dta_ReliefSurroundingVoted2016VsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVoted2016VsEIP_long$did <- dta_ReliefSurroundingVoted2016VsEIP_long$Relief * dta_ReliefSurroundingVoted2016VsEIP_long$d18

table(dta_ReliefSurroundingVoted2016VsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVoted2016VsEIPSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted2016VsEIP_long)

# summary(m_ReliefVoted2016VsEIPSurrounding)

m_ReliefVoted2016VsEIPSurrounding <- coeftest(m_ReliefVoted2016VsEIPSurrounding, vcov=vcovHC(m_ReliefVoted2016VsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVoted2016VsEIP_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.90747 -0.34080  0.01898  0.30259  1.01898 
# 
# Coefficients:
#   Estimate Std. Error  t value Pr(>|t|)    
# (Intercept)  0.6836367  0.0014567  469.306  < 2e-16 ***
#   d18         -0.5067270  0.0011880 -426.539  < 2e-16 ***
#   Relief      -0.0479336  0.0012543  -38.216  < 2e-16 ***
#   did          0.0674272  0.0017647   38.210  < 2e-16 ***
#   AgeCat20182  0.0587951  0.0013516   43.499  < 2e-16 ***
#   AgeCat20183  0.1894142  0.0012805  147.924  < 2e-16 ***
#   AgeCat20184  0.2166800  0.0013859  156.346  < 2e-16 ***
#   gender2     -0.0156912  0.0008881  -17.668  < 2e-16 ***
#   Dem         -0.0450177  0.0010848  -41.497  < 2e-16 ***
#   NPA         -0.1459265  0.0013074 -111.614  < 2e-16 ***
#   Third       -0.0637463  0.0061816  -10.312  < 2e-16 ***
#   Black        0.0071576  0.0012628    5.668 1.45e-08 ***
#   Hispanic    -0.0061709  0.0027984   -2.205   0.0274 *  
#   Other       -0.0342739  0.0020680  -16.573  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4259 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.2745,	Adjusted R-squared:  0.2745 
# F-statistic: 2.736e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted (by any method) in 2016 were MORE likely to early-in-person in 2018 than were those living in counties not hit by michael who voted (by any method) in 2016


# Election Day 2016 vs. EIP Only 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingEDOnly2016VsEIP_long <- reshape(dta_ReliefSurrounding,
                                                     varying = c("EDOnly2016", "EIPOnly2018"),
                                                     v.names = "voted",
                                                     timevar = "year",
                                                     times = c("EDOnly2016", "EIPOnly2018"),
                                                     new.row.names = 1:946286,    #multiply the number of obs by 2
                                                     direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingEDOnly2016VsEIP_long$year2018 <- ifelse(dta_ReliefSurroundingEDOnly2016VsEIP_long$year=='EIPOnly2018', 2018, 2016)
dta_ReliefSurroundingEDOnly2016VsEIP_long$d18 <- ifelse(dta_ReliefSurroundingEDOnly2016VsEIP_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingEDOnly2016VsEIP_long <- data.table(dta_ReliefSurroundingEDOnly2016VsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingEDOnly2016VsEIP_long$did <- dta_ReliefSurroundingEDOnly2016VsEIP_long$Relief * dta_ReliefSurroundingEDOnly2016VsEIP_long$d18

table(dta_ReliefSurroundingEDOnly2016VsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefEDOnly2016VsEIPSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingEDOnly2016VsEIP_long)

# summary(m_ReliefEDOnly2016VsEIPSurrounding)

m_ReliefEDOnly2016VsEIPSurrounding <- coeftest(m_ReliefEDOnly2016VsEIPSurrounding, vcov=vcovHC(m_ReliefEDOnly2016VsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingEDOnly2016VsEIP_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.3568 -0.2845 -0.2361  0.6537  0.8900 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.2552468  0.0014941 170.836  < 2e-16 ***
#   d18         -0.0368682  0.0012185 -30.257  < 2e-16 ***
#   Relief      -0.0441624  0.0012865 -34.328  < 2e-16 ***
#   did          0.0786164  0.0018100  43.435  < 2e-16 ***
#   AgeCat20182  0.0325878  0.0013863  23.506  < 2e-16 ***
#   AgeCat20183  0.0918677  0.0013134  69.949  < 2e-16 ***
#   AgeCat20184  0.0556894  0.0014215  39.177  < 2e-16 ***
#   gender2      0.0015942  0.0009109   1.750 0.080092 .  
# Dem         -0.0264095  0.0011127 -23.735  < 2e-16 ***
#   NPA         -0.0759330  0.0013410 -56.625  < 2e-16 ***
#   Third       -0.0341009  0.0063403  -5.378 7.51e-08 ***
#   Black        0.0080992  0.0012952   6.253 4.02e-10 ***
#   Hispanic    -0.0095858  0.0028703  -3.340 0.000839 ***
#   Other       -0.0251378  0.0021211 -11.851  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4368 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.01312,	Adjusted R-squared:  0.01311 
# F-statistic: 961.3 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted on election day in 2016 were MORE likely to vote early-in-person in 2018 than were those living in counties not hit by michael who voted on election day in 2016


# VBM Only 2016 vs. Voted 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVBMVsVoted_long <- reshape(dta_ReliefSurrounding,
                                                varying = c("VBMOnly2016", "Voted2018"),
                                                v.names = "voted",
                                                timevar = "year",
                                                times = c("VBMOnly2016", "Voted2018"),
                                                new.row.names = 1:946286,    #multiply the number of obs by 2
                                                direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVBMVsVoted_long$year2018 <- ifelse(dta_ReliefSurroundingVBMVsVoted_long$year=='Voted2018', 2018, 2016)
dta_ReliefSurroundingVBMVsVoted_long$d18 <- ifelse(dta_ReliefSurroundingVBMVsVoted_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVBMVsVoted_long <- data.table(dta_ReliefSurroundingVBMVsVoted_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVBMVsVoted_long$did <- dta_ReliefSurroundingVBMVsVoted_long$Relief * dta_ReliefSurroundingVBMVsVoted_long$d18

#table(dta_ReliefSurroundingVBMVsVoted_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVBMVsVotedSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVBMVsVoted_long)

# summary(m_ReliefVBMVsVotedSurrounding)

m_ReliefVBMVsVotedSurrounding <- coeftest(m_ReliefVBMVsVotedSurrounding, vcov=vcovHC(m_ReliefVBMVsVotedSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVBMVsVoted_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.83921 -0.29004 -0.06949  0.32171  1.13010 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.0642123  0.0014098  45.548  < 2e-16 ***
#   d18          0.4841377  0.0011497 421.089  < 2e-16 ***
#   Relief      -0.0308299  0.0012139 -25.398  < 2e-16 ***
#   did         -0.0894118  0.0017078 -52.354  < 2e-16 ***
#   AgeCat20182  0.0644657  0.0013081  49.282  < 2e-16 ***
#   AgeCat20183  0.1757983  0.0012392 141.860  < 2e-16 ***
#   AgeCat20184  0.2716895  0.0013413 202.563  < 2e-16 ***
#   gender2     -0.0175035  0.0008595 -20.365  < 2e-16 ***
#   Dem         -0.0283590  0.0010499 -27.011  < 2e-16 ***
#   NPA         -0.1230074  0.0012653 -97.216  < 2e-16 ***
#   Third       -0.0421348  0.0059824  -7.043 1.88e-12 ***
#   Black        0.0191739  0.0012221  15.689  < 2e-16 ***
#   Hispanic    -0.0190555  0.0027083  -7.036 1.98e-12 ***
#   Other       -0.0229688  0.0020014 -11.476  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4122 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.2718,	Adjusted R-squared:  0.2718 
# F-statistic: 2.699e+04 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted by mail in 2016 were less likely to vote (by any method) in 2018 than were those living in counties not hit by michael who voted by mail in 2016


# VBM Only 2016 vs. VBM Only 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVBM2016VsVBM2018_long <- reshape(dta_ReliefSurrounding,
                                                      varying = c("VBMOnly2016", "VBMOnly2018"),
                                                      v.names = "voted",
                                                      timevar = "year",
                                                      times = c("VBMOnly2016", "VBMOnly2018"),
                                                      new.row.names = 1:946286,    #multiply the number of obs by 2
                                                      direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVBM2016VsVBM2018_long$year2018 <- ifelse(dta_ReliefSurroundingVBM2016VsVBM2018_long$year=='VBMOnly2018', 2018, 2016)
dta_ReliefSurroundingVBM2016VsVBM2018_long$d18 <- ifelse(dta_ReliefSurroundingVBM2016VsVBM2018_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVBM2016VsVBM2018_long <- data.table(dta_ReliefSurroundingVBM2016VsVBM2018_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVBM2016VsVBM2018_long$did <- dta_ReliefSurroundingVBM2016VsVBM2018_long$Relief * dta_ReliefSurroundingVBM2016VsVBM2018_long$d18

table(dta_ReliefSurroundingVBM2016VsVBM2018_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVBM2016VsVBM2018Surrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVBM2016VsVBM2018_long)

# summary(m_ReliefVBM2016VsVBM2018Surrounding)

m_ReliefVBM2016VsVBM2018Surrounding <- coeftest(m_ReliefVBM2016VsVBM2018Surrounding, vcov=vcovHC(m_ReliefVBM2016VsVBM2018Surrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVBM2016VsVBM2018_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.29267 -0.14745 -0.10070 -0.06078  1.00752 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.1082075  0.0011434  94.638  < 2e-16 ***
#   d18         -0.0178466  0.0009325 -19.139  < 2e-16 ***
#   Relief      -0.0188105  0.0009845 -19.106  < 2e-16 ***
#   did         -0.0186397  0.0013851 -13.457  < 2e-16 ***
#   AgeCat20182  0.0162704  0.0010609  15.336  < 2e-16 ***
#   AgeCat20183  0.0580571  0.0010051  57.764  < 2e-16 ***
#   AgeCat20184  0.1839781  0.0010878 169.127  < 2e-16 ***
#   gender2     -0.0185291  0.0006971 -26.581  < 2e-16 ***
#   Dem         -0.0049673  0.0008515  -5.834 5.43e-09 ***
#   NPA         -0.0377343  0.0010262 -36.771  < 2e-16 ***
#   Third       -0.0064438  0.0048520  -1.328 0.184151    
# Black       -0.0034370  0.0009912  -3.468 0.000525 ***
#   Hispanic    -0.0041676  0.0021965  -1.897 0.057778 .  
# Other        0.0004833  0.0016232   0.298 0.765907    
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3343 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.04755,	Adjusted R-squared:  0.04753 
# F-statistic:  3610 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted by mail in 2016 were less likely to vote by mail in 2018 than were those living in counties not hit by michael who voted by mail in 2016


# VBM Only 2016 vs. EIP Only 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingVBMVsEIP_long <- reshape(dta_ReliefSurrounding,
                                              varying = c("VBMOnly2016", "EIPOnly2018"),
                                              v.names = "voted",
                                              timevar = "year",
                                              times = c("VBMOnly2016", "EIPOnly2018"),
                                              new.row.names = 1:946286,    #multiply the number of obs by 2
                                              direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingVBMVsEIP_long$year2018 <- ifelse(dta_ReliefSurroundingVBMVsEIP_long$year=='EIPOnly2018', 2018, 2016)
dta_ReliefSurroundingVBMVsEIP_long$d18 <- ifelse(dta_ReliefSurroundingVBMVsEIP_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingVBMVsEIP_long <- data.table(dta_ReliefSurroundingVBMVsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingVBMVsEIP_long$did <- dta_ReliefSurroundingVBMVsEIP_long$Relief * dta_ReliefSurroundingVBMVsEIP_long$d18

# table(dta_ReliefSurroundingVBMVsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefVBMVsEIPSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVBMVsEIP_long)

# summary(m_ReliefVBMVsEIPSurrounding)

m_ReliefVBMVsEIPSurrounding <- coeftest(m_ReliefVBMVsEIPSurrounding, vcov=vcovHC(m_ReliefVBMVsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingVBMVsEIP_long)
# 
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.41301 -0.23972 -0.17069 -0.05028  1.01387 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.0986793  0.0013427  73.493  < 2e-16 ***
#   d18          0.0896726  0.0010950  81.891  < 2e-16 ***
#   Relief      -0.0220523  0.0011561 -19.074  < 2e-16 ***
#   did          0.0502690  0.0016266  30.905  < 2e-16 ***
#   AgeCat20182  0.0298631  0.0012459  23.970  < 2e-16 ***
#   AgeCat20183  0.1091343  0.0011803  92.465  < 2e-16 ***
#   AgeCat20184  0.1784731  0.0012774 139.711  < 2e-16 ***
#   gender2     -0.0096001  0.0008186 -11.727  < 2e-16 ***
#   Dem         -0.0260342  0.0009999 -26.036  < 2e-16 ***
#   NPA         -0.0741845  0.0012051 -61.559  < 2e-16 ***
#   Third       -0.0314731  0.0056978  -5.524 3.32e-08 ***
#   Black        0.0179732  0.0011640  15.441  < 2e-16 ***
#   Hispanic    -0.0045110  0.0025794  -1.749 0.080319 .  
# Other       -0.0067081  0.0019062  -3.519 0.000433 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.3926 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.0549,	Adjusted R-squared:  0.05489 
# F-statistic:  4201 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted by mail in 2016 were MORE likely to vote early-in-person in 2018 than were those living in counties not hit by michael who voted by mail in 2016


# EIP Only 2016 vs. EIP Only 2018

# Reshape dataset from wide to long
dta_ReliefSurroundingEIPVsEIP_long <- reshape(dta_ReliefSurrounding,
                                              varying = c("EIPOnly2016", "EIPOnly2018"),
                                              v.names = "voted",
                                              timevar = "year",
                                              times = c("EIPOnly2016", "EIPOnly2018"),
                                              new.row.names = 1:946286,    #multiply the number of obs by 2
                                              direction = "long")



# Create year variable: 2018 = 1, 2016 = 0 
dta_ReliefSurroundingEIPVsEIP_long$year2018 <- ifelse(dta_ReliefSurroundingEIPVsEIP_long$year=='EIPOnly2018', 2018, 2016)
dta_ReliefSurroundingEIPVsEIP_long$d18 <- ifelse(dta_ReliefSurroundingEIPVsEIP_long$year2018==2018, 1, 0)

# Convert dataset to data.table
dta_ReliefSurroundingEIPVsEIP_long <- data.table(dta_ReliefSurroundingEIPVsEIP_long)


# DiD variable: treatment*year
dta_ReliefSurroundingEIPVsEIP_long$did <- dta_ReliefSurroundingEIPVsEIP_long$Relief * dta_ReliefSurroundingEIPVsEIP_long$d18

# table(dta_ReliefSurroundingEIPVsEIP_long$did)



# Difference-in-Difference County Treatment Model Estimation 
m_ReliefEIPVsEIPSurrounding <- lm(voted ~ d18 + Relief + did + AgeCat2018 + gender2 + Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingEIPVsEIP_long)

# summary(m_ReliefEIPVsEIPSurrounding)

m_ReliefEIPVsEIPSurrounding <- coeftest(m_ReliefEIPVsEIPSurrounding, vcov=vcovHC(m_ReliefEIPVsEIPSurrounding, type="HC0", cluster="CountyCode.x")) 

# Call:
#   lm(formula = voted ~ d18 + Relief + did + AgeCat2018 + gender2 + 
#        Dem + NPA + Third + Black + Hispanic + Other, data = dta_ReliefSurroundingEIPVsEIP_long)
# 
# Residuals:
#   Min      1Q  Median      3Q     Max 
# -0.4379 -0.3322 -0.2242  0.5918  0.9436 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.2620006  0.0015282 171.441  < 2e-16 ***
#   d18         -0.0771095  0.0012463 -61.869  < 2e-16 ***
#   Relief      -0.0010990  0.0013159  -0.835    0.404    
# did          0.0274050  0.0018513  14.803  < 2e-16 ***
#   AgeCat20182  0.0392869  0.0014180  27.706  < 2e-16 ***
#   AgeCat20183  0.1522662  0.0013434 113.347  < 2e-16 ***
#   AgeCat20184  0.1555562  0.0014540 106.988  < 2e-16 ***
#   gender2     -0.0049238  0.0009317  -5.285 1.26e-07 ***
#   Dem         -0.0340373  0.0011381 -29.907  < 2e-16 ***
#   NPA         -0.1044403  0.0013716 -76.144  < 2e-16 ***
#   Third       -0.0505844  0.0064851  -7.800 6.19e-15 ***
#   Black        0.0202902  0.0013248  15.315  < 2e-16 ***
#   Hispanic     0.0032636  0.0029358   1.112    0.266    
# Other       -0.0190862  0.0021696  -8.797  < 2e-16 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 0.4468 on 940112 degrees of freedom
# (6160 observations deleted due to missingness)
# Multiple R-squared:  0.03768,	Adjusted R-squared:  0.03767 
# F-statistic:  2832 on 13 and 940112 DF,  p-value: < 2.2e-16

# those living in counties hit by michael who voted early-in-person in 2016 were MORE likely to vote early-in-person in 2018 than were those living in counties not hit by michael who voted early-in-person in 2016























# Specify the width of your confidence intervals
interval1 <- -qnorm((1-0.9)/2)   # 90% CI
interval2 <- -qnorm((1-0.95)/2)  # 95% CI


# Put model estimates into temporary data.frames:

#================#
# Relief Policy  #
#================#

m_ReliefSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                         Coefficient = summary(m_ReliefSurrounding)$coef[4, 1],
                                         SE = summary(m_ReliefSurrounding)$coef[4, 2],
                                         modelName = "Voted 2016 vs. Voted 2018")

m_ReliefVBMVsVotedSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                   Coefficient = summary(m_ReliefVBMVsVotedSurrounding)$coef[4, 1],
                                                   SE = summary(m_ReliefVBMVsVotedSurrounding)$coef[4, 2],
                                                   modelName = "VBM 2016 vs. Voted 2018")

m_ReliefEDOnly2016VsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                        Coefficient = summary(m_ReliefEDOnly2016VsEIPSurrounding)$coef[4, 1],
                                                        SE = summary(m_ReliefEDOnly2016VsEIPSurrounding)$coef[4, 2],
                                                        modelName = "ED 2016 vs. EIP 2018")

# m_ReliefAnyEarlyVsVotedSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
#                                                       Coefficient = summary(m_ReliefAnyEarlyVsVotedSurrounding)$coef[4, 1],
#                                                       SE = summary(m_ReliefAnyEarlyVsVotedSurrounding)$coef[4, 2],
#                                                       modelName = "Any Early 2016 vs. Voted 2018")

m_ReliefEIPVsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                 Coefficient = summary(m_ReliefEIPVsEIPSurrounding)$coef[4, 1],
                                                 SE = summary(m_ReliefEIPVsEIPSurrounding)$coef[4, 2],
                                                 modelName = "EIP 2016 vs. EIP 2018")

m_ReliefVBM2016VsVBM2018Surrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                         Coefficient = summary(m_ReliefVBM2016VsVBM2018Surrounding)$coef[4, 1],
                                                         SE = summary(m_ReliefVBM2016VsVBM2018Surrounding)$coef[4, 2],
                                                         modelName = "VBM 2016 vs. VBM 2018")

m_ReliefVBMVsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                 Coefficient = summary(m_ReliefVBMVsEIPSurrounding)$coef[4, 1],
                                                 SE = summary(m_ReliefVBMVsEIPSurrounding)$coef[4, 2],
                                                 modelName = "VBM 2016 vs. EIP 2018")

m_ReliefVoted16VsAnyEarlySurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                          Coefficient = summary(m_ReliefVoted16VsAnyEarlySurrounding)$coef[4, 1],
                                                          SE = summary(m_ReliefVoted16VsAnyEarlySurrounding)$coef[4, 2],
                                                          modelName = "Voted 2016 vs. Any Early 2018")

m_ReliefVoted16VsEDSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                    Coefficient = summary(m_ReliefVoted16VsEDSurrounding)$coef[4, 1],
                                                    SE = summary(m_ReliefVoted16VsEDSurrounding)$coef[4, 2],
                                                    modelName = "Voted 2016 vs. ED 2018")

m_ReliefVoted16VsVBMSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                     Coefficient = summary(m_ReliefVoted16VsVBMSurrounding)$coef[4, 1],
                                                     SE = summary(m_ReliefVoted16VsVBMSurrounding)$coef[4, 2],
                                                     modelName = "Voted 2016 vs. VBM 2018")

m_ReliefVoted16VsEIPSurrounding_frame1 <- data.frame(Variable = "Differences-in-Differences",
                                                     Coefficient = summary(m_ReliefVoted2016VsEIPSurrounding)$coef[4, 1],
                                                     SE = summary(m_ReliefVoted2016VsEIPSurrounding)$coef[4, 2],
                                                     modelName = "Voted 2016 vs. EIP 2018")

allModelFrame_ReliefSurrounding <- data.frame(rbind(m_ReliefEDOnly2016VsEIPSurrounding_frame1, m_ReliefEIPVsEIPSurrounding_frame1,
                                                    m_ReliefVBMVsEIPSurrounding_frame1, m_ReliefVBM2016VsVBM2018Surrounding_frame1,
                                                    m_ReliefVBMVsVotedSurrounding_frame1, m_ReliefVoted16VsAnyEarlySurrounding_frame1,
                                                    m_ReliefVoted16VsEDSurrounding_frame1, m_ReliefVoted16VsEIPSurrounding_frame1,
                                                    m_ReliefVoted16VsVBMSurrounding_frame1, m_ReliefSurrounding_frame1))

# Plot
p_ReliefSurrounding <- ggplot(allModelFrame_ReliefSurrounding, aes(colour = Variable))
p_ReliefSurrounding <- p_ReliefSurrounding + geom_hline(yintercept = 0, colour = gray(1/2), lty = 2)
p_ReliefSurrounding <- p_ReliefSurrounding + geom_linerange(aes(x = Variable, ymin = Coefficient - SE*interval1,
                                                                ymax = Coefficient + SE*interval1),
                                                            lwd = 1, position = position_dodge(width = 1))

p_ReliefSurrounding <- p_ReliefSurrounding + geom_pointrange(aes(x = Variable, y = Coefficient, ymin = Coefficient - SE*interval2,
                                                                 ymax = Coefficient + SE*interval2),
                                                             lwd = 0.5, position = position_dodge(width = 1),
                                                             shape = allModelFrame_ReliefSurrounding$Variable)
#p_ReliefSurrounding <- p_ReliefSurrounding + coord_flip() + theme_bw()
# p_ReliefSurrounding <- p_ReliefSurrounding + facet_grid(modelName ~. , rows = 4, cols = 2)
p_ReliefSurrounding <- p_ReliefSurrounding + facet_wrap(~modelName, nrow = 4)
p_ReliefSurrounding <- p_ReliefSurrounding + theme_bw(base_line_size = 0.2) + ggtitle("") + theme(plot.title = element_text(hjust = 0.5))
p_ReliefSurrounding <- p_ReliefSurrounding + theme(legend.title=element_blank(), 
                                                   axis.title.x=element_blank(),
                                                   axis.text.x=element_blank(),
                                                   axis.ticks.x=element_blank()) 
p_ReliefSurrounding <- p_ReliefSurrounding + scale_y_continuous() + ylab("")
p_ReliefSurrounding <- p_ReliefSurrounding + scale_color_manual(values=c('#000000','#000000', '#000000', '#000000')) 
#p_ReliefSurrounding <- p_ReliefSurrounding + guides(colour = guide_legend(override.aes = list(shap_ReliefSurroundinge = c(1,2,3,4)),
   #                                                                       ncol = 2, nrow = 2))
p_ReliefSurrounding <- p_ReliefSurrounding + theme(legend.position='bottom', legend.box.background = element_rect(colour = "black"))
p_ReliefSurrounding <- p_ReliefSurrounding + theme(text = element_text(size=20))
#print(p_ReliefSurrounding)

ggsave("Relief.png", height=10, width=15, units='in', dpi=1000)






#=====================#
# Model Output/ LaTeX #  
#=====================#

library(lmtest) ; library(sandwich)
library(texreg)
library(xtable)


# Full DiD Models for Relief
ReliefModelsPart1 <- 
  texreg(list(m_ReliefSurrounding, m_ReliefEDOnly2016VsEIPSurrounding, m_ReliefEIPVsEIPSurrounding, m_ReliefVBM2016VsVBM2018Surrounding, m_ReliefVBMVsEIPSurrounding),
         caption="Relief Full DiD 2016-2018 Part 1",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("Voted 2016 Vs Voted 2018", "ED 2016 Vs EIP 2018", "EIP 2016 Vs EIP 2018", "VBM 2016 Vs VBM 2018", "VBM 2016 Vs EIP 2018"),
         override.se=list(summary(m_ReliefSurrounding)$coef[,2],
                          summary(m_ReliefEDOnly2016VsEIPSurrounding)$coef[,2],
                          summary(m_ReliefEIPVsEIPSurrounding)$coef[,2],
                          summary(m_ReliefVBM2016VsVBM2018Surrounding)$coef[,2],
                          summary(m_ReliefVBMVsEIPSurrounding)$coef[,2],
                          override.pval=list(summary(m_ReliefVBMVsEIPSurrounding)$coef[,4],
                                             summary(m_ReliefSurrounding)$coef[,4],
                                             summary(m_ReliefEDOnly2016VsEIPSurrounding)$coef[,4],
                                             summary(m_ReliefEIPVsEIPSurrounding)$coef[,4],
                                             summary(m_ReliefVBM2016VsVBM2018Surrounding)$coef[,4])))

ReliefModelsPart2 <- 
  texreg(list(m_ReliefVBMVsVotedSurrounding, m_ReliefVoted16VsAnyEarlySurrounding, m_ReliefVoted16VsEDSurrounding, m_ReliefVoted16VsVBMSurrounding, m_ReliefVoted2016VsEIPSurrounding),
         caption="Relief Full DiD 2016-2018 Part 2",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("VBM 2016 Vs Voted 2018", "Voted 2016 Vs Any Early 2018", "Voted 2016 Vs ED 2018", "Voted 2016 Vs VBM 2018", "Voted 2016 Vs EIP 2018"),
         override.se=list(summary(m_ReliefVBMVsVotedSurrounding)$coef[,2],
                          summary(m_ReliefVoted16VsAnyEarlySurrounding)$coef[,2],
                          summary(m_ReliefVoted16VsEDSurrounding)$coef[,2],
                          summary(m_ReliefVoted16VsVBMSurrounding)$coef[,2],
                          summary(m_ReliefVoted2016VsEIPSurrounding)$coef[,2],
                          override.pval=list(summary(m_ReliefVBMVsVotedSurrounding)$coef[,4],
                                             summary(m_ReliefVoted16VsAnyEarlySurrounding)$coef[,4],
                                             summary(m_ReliefVoted16VsEDSurrounding)$coef[,4],
                                             summary(m_ReliefVoted16VsVBMSurrounding)$coef[,4],
                                             summary(m_ReliefVoted2016VsEIPSurrounding)$coef[,4])))

# Full DiD Models for Relief (Robust Standard Errors)

ReliefModelsPart1 <- 
  texreg(list(m_ReliefSurrounding, m_ReliefEDOnly2016VsEIPSurrounding, m_ReliefEIPVsEIPSurrounding, m_ReliefVBM2016VsVBM2018Surrounding, m_ReliefVBMVsEIPSurrounding),
         caption="Relief Full DiD 2016-2018 Part 1",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("Voted 2016 Vs Voted 2018", "ED 2016 Vs EIP 2018", "EIP 2016 Vs EIP 2018", "VBM 2016 Vs VBM 2018", "VBM 2016 Vs EIP 2018"),
         override.se=list(m_ReliefSurrounding[,2],
                          m_ReliefEDOnly2016VsEIPSurrounding[,2],
                          m_ReliefEIPVsEIPSurrounding[,2],
                          m_ReliefVBM2016VsVBM2018Surrounding[,2],
                          m_ReliefVBMVsEIPSurrounding[,2],
                          override.pval=list(m_ReliefVBMVsEIPSurrounding[,4],
                                             m_ReliefSurrounding[,4],
                                             m_ReliefEDOnly2016VsEIPSurrounding[,4],
                                             m_ReliefEIPVsEIPSurrounding[,4],
                                             m_ReliefVBM2016VsVBM2018Surrounding[,4])))

ReliefModelsPart2 <- 
  texreg(list(m_ReliefVBMVsVotedSurrounding, m_ReliefVoted16VsAnyEarlySurrounding, m_ReliefVoted16VsEDSurrounding, m_ReliefVoted16VsVBMSurrounding, m_ReliefVoted2016VsEIPSurrounding),
         caption="Relief Full DiD 2016-2018 Part 2",
         digits = 3,
         dcolumn=FALSE,
         model.names=c("VBM 2016 Vs Voted 2018", "Voted 2016 Vs Any Early 2018", "Voted 2016 Vs ED 2018", "Voted 2016 Vs VBM 2018", "Voted 2016 Vs EIP 2018"),
         override.se=list(m_ReliefVBMVsVotedSurrounding[,2],
                          m_ReliefVoted16VsAnyEarlySurrounding[,2],
                          m_ReliefVoted16VsEDSurrounding[,2],
                          m_ReliefVoted16VsVBMSurrounding[,2],
                          m_ReliefVoted2016VsEIPSurrounding[,2],
                          override.pval=list(m_ReliefVBMVsVotedSurrounding[,4],
                                             m_ReliefVoted16VsAnyEarlySurrounding[,4],
                                             m_ReliefVoted16VsEDSurrounding[,4],
                                             m_ReliefVoted16VsVBMSurrounding[,4],
                                             m_ReliefVoted2016VsEIPSurrounding[,4])))

